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1.  Introduction  1 I I 

Let  X 3 [X^:t  >0}  be  a regenerative  process  which  we  wish  to 
simulate.  Under  mild  regularity  conditions  the  distribution  of  con- 

verges to  the  distribution  of  some  limiting  random  variable  (or  vector)  X. 
This  type  of  convergence  is  known  as  weak  convergence  and  written 
X^  ^ X,  as  t t 00.  Simulators  speak  of  X as  the  "steady- state”  con- 
figuration of  the  system  and  are  often  interested  in  estimating  the 
constant  r s E[f(X)},  where  f is  a given  real-valued  function  defined 
on  the  state- space  of  the  process  X.  The  regenerative  method  of  simula- 
tion provides  a means  of  constructing  point  and  interval  estimates  for  r; 
see  IGLEHART  (1977)  for  an  expository  summary  of  this  method. 

The  problem  we  consider  in  this  paper  does  not  involve  estimation 
of  r^  but  rather  the  estimation  of  extreme  values  of  the  regenerative 
process  X.  Suppose^  for  the  sake  of  discussion,  we  are  simulating  a 

stable  GI/G/1  queue  in  order  to  estimate  the  maximum  waiting  time  among 

* 

the  first  n+l  customers;  call  this  random  variable  W^.  As  n grows, 

so  will  W . However.  W does  not  converge  to  a finite  limit,  but  rather 
n ' n ' 
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diverges  to  -hh.  We  will  be  interested  In  estimating  the  distribution 
* 

function  of  W^  for  finite,  but  large  n.  By  the  same  token  we  might 
wish  to  estimate  the  maximum  queue  length  during  the  Interval  [0,t]. 

While  this  problem  of  estimating  extreme  values  would  seem  to  be  of  great 
practical  importance  to  simulators,  we  know  of  no  papers  In  the  simulation 
literature  which  offer  any  guidance  on  the  subject.  This  paper  will 
attempt  to  partially  fill  the  gap. 

We  begin  in  Section  2 by  summarizing  a series  of  probabilistic 
results  in  extreme  value  theory  which  will  provide  the  theoretical  basis 
for  the  methods  we  propose.  Section  3 discusses  two  methods  for  estimating 
extreme  values  for  the  general  regenerative  simulation.  In  Sections  k and  3 
we  treat  the  special  cases  of  the  GI/G/1  queue  and  birth-death  processes 
respectively.  Theoretical  results  are  available  for  these  two  classes 
of  regenerative  processes  that  are  useful  in  assessing  the  accuracy  of 
the  simulation  methods  proposed.  Section  6 contains  the  numerical  results 
for  simulations  of  the  a/H/ I queue  carried  out  to  illustrate  the  estima- 
tion methods  proposed. 


2.  Probabilistic  Background 

Let  (F^:n  > 1)  be  a sequence  of  distribution  functions  (d.f. 's) 

on  the  real  line,  IR  = (-»,  +oo) . This  sequence  converges  weakly  to  a 

d.f.  G if  lim  which  are  continuity 

n -♦  00 

points  of  G.  We  write  ^ G to  denote  this  type  of  convergence.  If 
(resp.  X)  is  a random  variable  (r.v.)  with  d.f.  F^  (resp.  G) , we 
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also  write  X to  denote  this  weak  convergence.  Sometimes  It  is 

n 

convenient  to  write  X^  ^ G to  connote  the  same  thing.  The  material 
presented  In  this  section  can  be  found  for  the  most  part  In  deHAAN  (1970), 
currently  the  best  comprehensive  treatment  of  Che  subject. 

Now  let  ^ be  a sequence  of  Independent,  Identically 

i distributed  (l.l.d.)  r.v. 's  and  denote  Che  maximum  of  the  first  n r.v.'s 

by  3 max(X^: 1 

will  have  d.f.  We  shall  say  chat  F belongs  Co  Che  domain  of 

attraction  of  Che  non-degenerate  d.f.  G,  and  write  F £ i)(G),  If  we  can 

^ choose  two  sequences  of  constants  ^ with 

a > 0 such  that 
n 


< k < n).  If  each  of  Che  then 


H ! 

hk 


(2.1) 


f 'V  * ■’n) 


as  n -»oa  for  all  x €11  for  which  G Is  continuous.  Equivalently, 


F € fl(G)  if  n -» <».  Thus  for  large  n we  would 


approximate  P{M^  < x)  by  G( ( x-b^)/a^) , If  a r.v.  X has  d.f. 


F € i)(G),  we  also  write  X € i)(G). 

A famous  result  In  extreme  value  theory  states  that  Che  only  d.f. 's 
G which  can  arise  In  (2.1)  are  of  one  of  Che  following  three  types: 


X < 0 


(2.2) 


exp(  -x”^  , X > 0 


(2.3) 


= 


exp{.(.x)  ) , X < 0 


(2.i^) 


A(x)  - exp(-e'*)  , 


X > 0 


X en  , 


where  in  (2.2)  and  (2.3)  OC  Is  a poslclve  constant.  Recall  that  two 
d.f. 's  and  are  said  to  be  of  the  same  type  if  there  exists 

two  constant  a and  b,  a > 0,  such  that  G^(x)  s C2(ax  + b)  for  all 
X €]R.  Thus  aside  from  translations  and  scaling  by  a positive  constant 
the  three  d.f. *3  given  In  (2.2)  - (2.4)  are  the  only  ones  that  can  appear 
In  (2.1).  This  result  on  the  three  types  of  limit  d.f.'s  Is  usually 
attributed  to  GNEDENKO  (19^3),  however  It  was  first  formulated  In  this 
way  by  FISHER  and  TIPPETT  ( I923) . 

The  next  logical  result  to  seek  Is  necessary  and  sufficient  conditions 

for  F £ D(G),  where  G of  necessity  Is  one  of  the  three  d.f.  's  given  In 

(2.2)  - (2.4).  Furthermore,  If  F € )0(G)  we  need  a method  for  selecting 

the  two  sequences  fa  :n  > 1)  and  (b  :n  > 1).  To  this  end  we  first 

n ~ n — 

define  the  right  endpoint,  x < +00,  of  the  d.  f . F as 


x^  = sup(x:F(x)  < 1) 


A d.f.  F £ only  If  for  all  x > 0 


(2.5) 


, , 1 - F(  tx)  -OC 

1 ■ F(t)'  = * 

t -»  00  ' ' 


I 


If  F e then  we  can  Cake 


(2.6)  » inf(x:l  - F(x)  < l/n) 

and  = 0.  A d.f.  F € if  if  *o  ^ “ *“** 

X > 0 


(2.7) 


1 - FIXq  - (Cx)‘  ] 
lim  . ■ - X 


t -»  00 


i - F(t) 


If  F € , then  we  can  take  » x^  and 


a^  = Xq  - inf{x:l  - F(x)  < l/n) 

The  final  case^  F € il(A),  is  the  most  important  one  for  our 
applications.  A d. f . F € ^(A)  if  and  only  if 


(2.8) 


lim  1 - F(t  + xf(t))  -X 
t t Xq  1 - F(  c) 


where  for  t < x^ 


for  all 


f(t) 


(1  - F(s))ds 
1 - F(  C) 


If  F € fl(A),  then  we  can  take 


for  all 


simulation 


X €3R  , 
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= lnf{x;l  - F(x)  < 1/n) 

and 

*0 

/ (I  - F(c)Jdc 
*n  “ I - F(b  ) • 

Alternative  expressions  are  available  for  a and  b . Let  Q 

n n ^n' 

denote  the  p- quantile  of  the  d.f.  f":  for  0 < p < 1^ 

Qn(p)  = inf{x:F"(x)  > p)  . 

Then  if  F € i)(A),  we  can  alternatively  select 

(2.9)  = 

and 

(2.10)  . Q_^(e-*)  . 


Furthermore,  if  F € jO(A)  and  x.  = +«  M /a  =»  1 as  n Many 

u ' n n 

of  the  classical  d.f. 's  such  as  the  exponential,  ganma,  normal,  lognormal, 
logistic,  and  Cauchy  belong  to  fl(A). 

Suppose  F has  x^  = +»  and  possesses  an  exponential  tail: 

(2.11)  1 - F(x)  b exp(-ax)  , as  x -»oo  , 
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where  a and  b are  two  positive  constants.  Then  It  is  easy  to  check 

that  (2.8)  holds  and  F € fi(A).  Using  the  expressions  (2.9)  and  (2.10) 

It  can  be  shown  that  b and  a can  be  selected  as  follows: 

n n 


b^  ■ a"^  In(nb)  , 

and 

-1 

a a a 
n 


An  Interesting  (and  practical)  situation  arises  if  F is  a disctJte 
d.f.  as,  for  example,  the  geometric  d.f.  F(x)  » 1 - exp(-(x]),  x > 0, 
where  [x]  is  the  integer  part  of  x.  In  this  case  neither  (2.5)  nor 
(2.8)  hold,  and  since  = +oo,  F does  not  belong  to  the  domain  of 
attraction  of  any  of  the  three  types  (2.2)  - (2.I4.).  However,  a result 
has  been  salvaged  by  ANDERSON  ( I97O) . Let  a be  the  class  of  all  d.f. 's 
whose  support  consists  of  all  sufficiently  large  positive  integers.  Then 
if  F € a. 


(2.12) 

and 

(2.13) 


lim  sup  f''(0j"^x  + b ) < exp(-e”*) 
n -» 00 


lim  inf  f"(Q£"^x  + b ) > exp(-e”^*”^^ ) 
n -»  00 


for  some  a > 0,  all  x,  and  some  sequence  (b  :n  > 1]  i'f  and  only  if 


(2.1k) 


1 - Ff  n)  a 

1 . F(nUi  = ® 


n -♦  00 
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When  this  condition  holds.  Che  constants  b can  be  selected  as  follows. 

• n 


For  F 6 G and  each  positive  integer  n let  h(n)  ■ -log(l-F(n))  and 

define  h to  be  the  extension  of  h obtained  by  linear  interpolation 
c 

for  X > 1.  Then  define  for  x > 1 

F^(x)  = 1 - exp(-h^(x))  . 

Clearly  F^  is  a continuous  d. f . and  for  sufficiently  large  x is 
strictly  increasing  since  F £ G.  For  x < 1 the  F^  can  be  defined 
arbitrarily  just  so  long  as  it  is  a d.f.  In  terms  of  F^  we  can  define 
b^  for  large  n as  Che  unique  root  of 

I 

I 

I 

1 - F (b  ) = 1/n  . 

c'  n' 

If  F € G and  1 - F(n)  ~ b exp(-an)  as  n -»<»  (a,b  > 0)^  then 

for  b = a”^  In(nb)  it  can  easily  be  shown  using  Che  method  followed 
n 

by  HEYDE  ( I97O)  that  for  integer  k 

(2.15)  lim  - (bj  < k}  - exp( -e"®^ ) 1 _ 0 , 

n — » 00  L 4 

where  d = b - [b  1.  Thus  for  n large  we  would  use  Che  approximation 
n n '■  n^ 


I 

P(M^  < k + [bj  ) ar  exp(.e-^^‘‘-‘*n))  ^ 

P(M^  < k)  ? exp(.e'®(‘‘'‘’"^)  . 

Suppose  now  that  we  also  have  defined  on  Che  probability  triple 
(0,j,P)  Chat  supports  the  i.i.d.  sequence  (X^:n  > 1)  a renewal  process 
(i(t):t  > 0)  with  mean  time  between  renewals  m (0  < m < o»).  Then 
the  weak  law  of  large  numbers  for  renewal  processes  states  that 
i5(t)/  t ^ m ^ as  n -» oo.  Next  set 


or 

(2.16) 


= max(X^: 1 < k < £( t) } . 

The  following  useful  result  for  this  situation  was  obtained  by  BERMAN 
(1962).  If  ^ three  extreme  value  d.  f. 's 

(2.2)  - (2,k)j  Chen  as  t -» « 

(2.17)  • 

This  result  provides  a useful  tool  for  extreme  values  of  regenerative 
processes.  To  be  explicit  suppose  X = {X^:t  > Oj  is  a regenerative 
process  defined  on  (!i,3,P)  and  Tj^,  k > 1^  is  the  time  of  the  kth 
regeneration  point  of  X with  T^  = 0.  Then  the  renewal  process 
(i(c):t  > 0)  which  counts  the  number  of  regeneration  points  in  (O^t] 
is  defined  by 
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i(  t)  3 max{k:Tj^  < c) 
with  £(0)  3 0.  For  k > let 

Since  X is  regenerative^  the  sequence  of  maxima^  **• 

i.i.d.  Then  if  = sup{X^;0  < s < t)^  clearly 

(2. 18)  max{Mj^:  1 < k < t)  ) < < raax{M^;  I < k < «( t)  ♦ 1)  . 

Combining  the  inequalities  of  (2.15)  with  the  limit  theorem  of  (2.1k) 
enables  us  to  show  that 

(2.19)  ^^t'^’ltj^/^t]  > 

where  m = provided  £ i)(G).  Of  course  if  € 0,  then  the 

weaker  results  of  Anderson  or  Heyde  are  all  that  can  be  expected. 

We  conclude  this  section  by  summarizing  the  problems  confronting 
us  for  a regenerative  processes  with  continuous  state  space.  If  € i)(G), 
then  we  can  use  (2.I9)  to  obtain  the  asymptotic  (for  large  t)  approxima- 
tion 


(2.20) 


P(Lj.  < x)  = G 


l/m 


(llAli) 


iLi 
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If  the  sinnilation  is  run  for  n cycles,  then  (2.20)  should  be  replaced 
by 


(2.21) 


I max  M^<xl3rG( 
' 1 < k < n ^ ^ 


For  (2.20)  or  (2.21)  to  be  useful,  we  must  estimate  m.  a , and  b . The 

expected  cycle  length,  m,  can  of  course  be  estimated  by  the  sample  mean 

of  the  cycle  lengths  observed.  Two  methods  for  estimating  a and  b 

n n 

will  be  discussed  in  Section  3*  Finally,  we  must  assess  whether 

£ i)(G)  for  one  of  the  three  d.f. 's  (G's)  given  in  (2.2)  - (2.U.). 

For  most  simulations  in  which  extreme  values  are  being  estimated,  the 

limit  d.f.  's  G will  be  either  A or  since  the  maxima  arising  are 

unbounded  {^q  = +”) • Our  experience  with  specific  examples  indicates 

that  if  the  regenerative  process  is  stable  (converges  to  a non-degenerats 

limit),  then  G = A.  While  if  the  process  is  "null-recurrent" 

(m  = E(Tj^}  = +*),  then  G = However,  in  this  case  G^™(x)  = 0 

for  all  x > 0 which  indicates  that  a different  normalization  must  be 

used  to  obtain  a non-degenerate  limit.  In  any  case,  we  note  that  if 

X £ '®(*^q)  with  constants  a^  > 0 and  b^  = 0,  then  In  X £ fl(A)  with 

constants  a*  = Of  ^ and  b ' = O:  ^ In  a , For  Che  balance  of  this  paper 
n n n 

we  shall  assume  that  G = A for  continuous  state  space  processes.  We 
note  in  passing  that  the  extreme  value  behavior  of  some  function  of  a 
regenerative  process  can  be  handled  in  Che  same  way.  If  Che  state  space 
of  the  regenerative  process  is  discrete,  then  we  shall  only  consider  Che 
situation  in  which  the  d.f.  of  £ Q and  P{M^  > n)  -x  b exp(-an) 
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as  n -» oo  for  some  a^b  > 0.  In  this  case  we  can  approximate  the  d.f. 
of  or  max{M^:l  < k < n)  by  using  (2.12)  and  (2.13)  or  (2.15) 

(2.16). 


. Statistical  Estimation  Problem 


We  have  seen  in  (2.20)  and  (2.21)  that  the  key  to  estimating  the 

d.f.  of  extreme  values  occurring  in  regenerative  processes  is  the  constants 

a and  b . These  constants  are  in  turn  determined  by  the  tail  behavior 
n n 

of  the  d.f.  of  the  maximum  value  of  the  process  within  a cycle. 

Following  Che  remarks  at  Che  end  of  Section  2^  we  shall  first  assume 

the  state  space  of  the  regenerative  process  is  continuous  and  that 

€ 5(A).  Two  methods  are  proposed  here  for  estimating  a^  and  b^. 

The  first  is  based  on  the  asymptotic  relation  (2.11). 

Assume  that  n cycles  are  simulated  and  the  n maxima^  .5  ^ 

are  arranged  in  decreasing  order.  Call  them  Y.  ; i.e.,  Y,  > Y_  > 

j,n  ^ ~ ~ 

•••  > Now  form  the  tail  of  the  empi  iCal  d.f.,  namely. 


E„(x) 


x)  = n i:  1 < i < n 


(i:l  < i < n,  Y.  > x} 
— — > in 


Plot  log  E^(x)  versus  x.  If  this  graph  is  roughly  linear  in  x (at 
least  for  large  x),  then  we  can  assume  (2.11)  holds  and  fit  a linear 
regression  line.  To  do  this  some  judgment  will  have  to  be  made  as  to 
how  large  x need  be  before  the  relationship  is  linear.  Using  the 
standard  point  and  interval  estimates  for  the  slope  and  y- intercept  of 


12 
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a regression  line  the  parameters  a and  b of  (2.11)  can  be  estimated 
by  a and  b^  say.  In  particular^  if  the  regression  y = cx-*-d  is 
fitted  to  the  plot  of  log  vs.  x,  then  a = -c  and  b * exp(d). 

This  in  turn  provides  estimates  for  In(nb)  and  a^  = a 

namely  b^  = In(nb)  and  Should  the  plot  of  log  E(x) 

vs.  X not  be  linear^  our  only  suggestion  is  to  first  take  logarithms 
of  the  and  try  again.  Maybe  the  underlying  situation  was 

m|  s cc  > 0. 

A second  method  for  estimating  a^  and  b^  is  that  given  by 

FEIGIN  and  IfElSSMAN  (1977)  which  is  based  on  the  k (a  fixed  positive 

integer)  largest  Y 's.  They  assume  that  m|  € f)(A).  Assuming  that 

n is  large  enough  so  that  ^ WEISSMAN  (197^) 

has  shown  that  the  asymptotic  maximum  likelihood  estimators  (AMLE)  of 

a and  b are  given  by 
n n 


a = Y,  - Y, 
n k.n  k.n  ^ 


/N  /S 

b = a In  k + Y, 
n n k.n 


and  the  asymptotic  uniformly  minimum  variance  unbiased  estimators 
(AUMVUE)  by 


a*  = Y,  , - Y, 

* * 
b = a I 

n k-l,n  k,n  > 

n n 

j 

k-1 

where  Y = (l/j)  Z Y 

isl 

j=l 

Euler's  constant.  Furthermore,  FEIGIN  and  WEISSMAN  (1977)  have  shown  that 


100(1-0:)^  (asymtocic)  confidence  intervals  for  a and  b are 

n n 

given  by 


and 

- I-  \n  - < . 

2 2 

where  p)  is  the  100p%  point  of  the  X distribution  with  r degrees 

of  freedom.  Also  is  the  distribution  of  the  ratio 

which  nij^  and  Y are  independent^  exp(-inj^)  has  a gamma  distribution 

with  parameters  k and  1^  and  (k-l)Y  is  gamma  with  parameters 

(k-1)  and  1.  The  quantiles^  U,  ,(p),  have  been  computed  by  FEIGIN 

and  WEISSMAN  (1977)  reproduced  in  Table  1. 

Assume  now  that  the  state  space  of  the  regenerative  proces  is 

discrete^  £ Q,  and  P(M^  > n)  — b exp(-an)  as  n -> « for  some  a^b  > 0. 

Then  we  can  use  the  approximation  given  in  (2.16),  ^ < k)  = A(a(k-b^)),. 

where  b^  = a In(bn).  Thus  again  we  can  apply  either  the  regression 

or  Feigin-Weissman  method  to  estimate  a and  b . 

n Cl 

Once  the  constants  a and  b have  been  estimated,  the  d.f.  of 

n n ' 

extreme  values  is  estimated  using  either  (2.20)  or  (2.21)  depending  on 
whether  the  simulation  was  run  for  a fixed  length  of  time,  t,  or  for  a 
fixed  number  of  cycles  n. 

These  methods,  which  seems  to  be  the  simplest,  for  estimating  a^ 
and  b^  will  be  illustrated  in  Section  6 for  the  M/M/1  queue.  Before 
leaving  this  section  we  point  out  some  other  relevant  references.  The 
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reliability  theory  literature  contains  many  references  on  the  problem 
of  testing  whether  observations  come  from  an  exponential  or  extreme 
value  d.f.  and  of  estimating  the  associated  parameters.  Two  useful 
places  to  find  such  papers  are  EPSTEIN  ( I960)  and  MANN,  SCHAFER,  and 
SINGPURWALLA  (I97U),  Chapter  5.  PICKANDS  (1975)  has  developed  a method 
for  determining  which  G d.f.  is  appropriate  for  a given  set  of  observa- 
tions. His  method  uses  a random  and  increasing  number  of  the  Y.  's 
as  n increases.  The  method  is  expensive  computationally  and  exphasizes 
an  aspect  of  the  extreme  value  problem  which  is  not  of  great  concern  for 
simulation. 


k.  GI/G/1  Queue 

The  GI/g/1  queue  and  the  birth-death  processes  treated  in  Section  5 

are  among  the  very  few  regenerative  processes  for  which  we  know  the 

values  of  a and  b . For  this  reason  these  processes  are  excellent 
n n 

candidates  for  testing  the  effectiveness  of  the  estimation  procedures 
proposed  in  Section  3. 

In  the  GI/g/ 1 queue  we  assume  customer  0 arrives  at  t^  = 0, 
finds  a free  server,  and  experiences  a service  time  v^.  Customer  n 
arrives  at  time  t and  experiences  a service  time  v . Customers  are 
served  in  their  order  of  arrival  and  the  server  is  never  idle  if  customers 
are  waiting.  Let  the  interarrival  times  t - t , = u n > 1.  We 
assume  the  two  sequences  ^ ®®ch  consist  of 

i.i.d.  r.v. 's  and  are  themselves  independent.  Let 

Efv  j = Li”^,  where  0 < ^.  u < ».  The  traffic  intensity  p = is 

n f 


!■  i 

Li 


[ assumed  to  be  less  than  one.  We  exclude  the  deterministic  system  in 

' which  both  the  v 's  and  u 's  are  degenerate.  Let  the  waiting  time  of 

n n 


the  nth  customer  be  W^^  the  workload  (or  virtual  waiting  time)  at 

time  t be  and  the  number  of  customers  in  the  system  at  time  t 

be  Q...  Also  set  W = max(W,  :0  < k < n],  V*  = sup(V  ;0  < s < t), 
t n k — — ^ t s — — ^ 

and  = sup(Q  :0  < s < t}.  Let  X = v ,-u,n>l.  and  set 
t *^’■^8  — — n n-1  n^  — ^ 

S = X,  + • • • + X for  n > 1 and  S.  =0,  If  n,  denotes  the  number 
n 1 n — 0 K 

of  customers  served  in  the  kth  busy  period^  then  n^  is  related  to  the 
partial  sum  process  since 

nj^  = inf{n  > 0:S^  < 0)  . 

When  p < 1^  we  have  m = E(n^}  < oo.  Also  -S  is  the  length  of  the 
first  idle  period.  We  assume  that  X^^  has  an  aperiodic  d.f.  (support 
is  not  concentrated  on  a set  of  points  of  the  form  0^  +h^  ±2^,  iJh, 
that  there  exists  a positive  number  k such  that  E{exp( (cX^^) ) = 1^ 
and  0 < p = E(X,  exp(KX,))  < «.  These  assumptions  will  normally  be 
satisfied  if  the  d.f.  of  v^  has  an  exponentially  decaying  tail;  e.g.^ 
when  v^  has  a gamma  distribution.  Under  these  conditions  we  know 
(see  IGLEHART  (1972))  that 

(It-.l)  (W*  - X ^ log  bj^n)/<"^  ^ A^™(x)  . 

and 

(1^.2)  (V*  - log  b2t)/x"^  A^"(x)  , 

where 


* » 

Thus  Co  use  (U-.l)  and  for  estimating  the  d.f. 's  of  W and 

we  need  only  estimate  m and  E{e  assuming  that  4^,  and 

/CVrt 

E(e  can  be  calculated  numerically.  In  Che  special  case  of  M/G/l 

- 1 kS 

queues  no  estimation  is  required^  since  m=  ( 1-p)  and  E(e  "1) 

= VCA+k).  If  the  simulation  is  carried  out  for  a fixed  number  of  cycles^ 
then  counterparts  of  (1<-.1)  and  (U.2)  hold  with  the  exponents  of  A 
removed. 

The  queue-length  process  > 0}  is  discrete-valued  and  the 

associated  d.f.  of  € 0.  Hence  a limit  theorem  comparable  to  (U-.l) 
or  (4.2)  does  not.  exist.  Instead  we  must  seek  results  like  (2.12)  and 
(2.13)  or  (2.15)  and  (2.l6).  Unfortunately^  these  results  are  only 
known  for  the  H/G/l  and  GI/M/1  queues,"  see  COHEN  (I968),  Theorems  7.2 
and  7.5.  Let  = sup{Q^:Tj^  ^ ^ ^ ^ M/G/l  queue 

the  counterpart  of  (2.16)  is 

f 

l 

-a(k-b  ) 

: (4.3)  PC  max  M < k)  = exp(-e  ) , 

[ 1 < j < " 

[ where 

i 

f 
! 

and 

1 I 

[■ 

i' 

(■ 

L . 


a = log((>»+<)/X) 
b^  = a‘^  log(b2n)  . 
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On  the  other  hand^  for  GI/M/ I queues  (4.5)  holds  with  a = log(  (^x-<)/u) 
and  the  same  value  for  b^.  Tables  2 and  5 contain  the  values  of 
bj^,  and  bg  for  the  M/M/ 1 and  M/E^/l  queues  as  a function  of 
the  traffic  intensity  p. 


5.  Birth-Death  Processes 

A second  class  of  regenerative  processes  for  which  theoretical 

results  are  available  is  birth-death  processes  in  discrete  or  continuous 

time.  Let  (X  :n  >0}  be  a discrete  time  Markov  chain  with  state-space 
n ^ 

E = (0^  1,  2^  . . . ) and  transition  probabilities  given  by 


(5.1) 


Pij  = Pi^ 


j = i-1 
j = i+1 
other 


where  ‘Iq  = Pq  “ ^ other  q^'s  and  p^'s  are  positive. 

This  chain  will  automatically  be  both  irreducible  and  periodic.  Further- 
more^ recall  that  it  will  be  recurrent  if  and  only  if 


Z (jr  p.)  = oo 

j=l  J 

where  JTq  ^ ^ j = ( Pq  ’ ' ’ Pj-l^^^'^l  ast,ume  the  chain 

is  recurrent.  It  will  be  positive  recurrent  if  and  only  if 
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p 

m 

< 

*>1 

•>2 

.1 

1.11 

9 

.900 

.09 

.9 

.2 

1.25 

8 

.1+00 

VO 

• 

.8 

.3 

l.i^3 

7 

CVJ. 

.21 

.7 

.k. 

1.67 

6 

.150 

.21+ 

.6 

.5 

2.00 

5 

.100 

.25 

.5 

.6 

2.50 

k 

.067 

.21+ 

.1+ 

.7 

3.33 

3 

.01+3 

.21 

.3 

.8 

5.00 

2 

.025 

.16 

.2 

.9 

10.00 

1 

.011 

.09 

.1 

• 95 

20.00 

0.5 

.005 

.01+75 

.05 

.99 

100.00 

0.1 

.001 

.0099 

.01 

TABLE  3 


Parameter  Values 

1 for 

M/Eg/l 

Queue  with 

V*  = 10 

p 

m 

K 

‘»2 

.1 

1.11 

15.00 

.3375 

.1562 

2.5 

.2 

1.25 

12.60 

.2016 

.2346 

1.7119 

1.^3 

10.61 

.1395 

.2874 

1.3038 

.k 

1.67 

8.83 

.1012 

.3179 

1.0201 

.5 

2.00 

7.19 

.0741 

.3263 

0.7957 

.6 

2.50 

5.64 

.053^^ 

.3118 

0.6050 

.7 

3.33 

4. 16 

.0367 

.^r33 

0.4357 

.S 

5.00 

2.73 

.022J 

.2094 

0.2809 

• 9 

10.00 

1.35 

.0106 

.1188 

0. 1366 

•95 

20.00 

0.67 

.0051 

.0630 

0.0674 

• 99 

100.00 

0.13 

.0010 

.0132 

0.0134 
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Next  define 


E JT  . < » 
j=o  J 


Tj^(k)  = inf(n  > 0:X^  = k). 


k € E , 


the  first  entrance  time  to  state  k.  Let  P^{*)  = P(*|Xq  = i}^  the 
conditional  probability  of  an  event^  given  = i.  Then  our  concern 
here  will  be  in  the  probability,  given  = i,  of  the  Markov  chain 

entering  state  n before  it  enters  state  0.  Let  this  probability  be 
denoted  by 


^iCn)  = P^{T^(n)  < T^(0)}  , 


i € (1,  2,  . . n-1]  . 


Fortunately,  this  probability  has  been  calculated  and  in  particular 

= (1  + S ; 

i=l 

see  CHUNG  (I96O),  p.  68.  Note  that  lim  = 0 when  the  chain  is 

n -*  00 

recurrent,  in  keeping  with  our  intuition.  Define 


= sup{X^:0  < n < Tq-1) 


Then 

(5.2)  Pq{M^  > n]  = rQ(n+l)  = ( ^ . 


3s  n "4  00  , 
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Suppose  now  that  we  have  a birth-death  process  (X^:t  > 0):  a 
continuous  time  Markov  chain  with  state  space  E = (0,  1^  ...)  and 
embedded  jump  chain  whose  probabilities  are  given  by  (5. 1).  As  above^ 
define  the  first  entrance  time  to  state  k and  the  maximum  in  the  first 
cycle  by 


-r^(k)  = inf(s  > 0:X^_  j,  X^  = j) 


and 


= sup{Xj.:  0 < t < Tj^(O) ) . 

Because  of  the  path  structure  of  the  birth-death  process,  it  is  easy  to 
see  from  (5-2)  that 

(5.5)  Po{M^  > n)  = [l  + , 

where  [resp.  are  the  birth  [resp,  death]  parameters  and 

”^0  ~ "^i  " ^^0^1  ”*  ^i-l^^l‘^2  '^i^  • argument  can  of 

course  be  used  to  show  that  (5.5)  also  holds  for  semi-Markov  processes 
with  embedded  jump  chain  whose  probabilities  are  given  by  (5.I). 

(5.4)  EXAMPLE.  I^M/s  queue.  The  queue- length  process,  > 0), 

is  a birth-death  process  with  parameters  = X and  = n (j  A s), 
j > 0.  Assume  the  queue  has  traffic  intensity  p = X/^s  < 1,  a necessary 
and  sufficient  condition  for  recurrence.  Then  from  (5-5) 


1 
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n .T I 


?oK  > »)  = [ j; 

U i=0  -I 


r s , , n-(s+l)  .-,-1 

^ 0 1 


S p isO 


Asymptotically^  as  n -+« 


it  > n)  ~ I 


(sVs.')n~^  , 


PqCm;  > n} 


I (s®(l- 


p)/s.')p  , 


p = 1 

p < 1 . 


Thus  for  p < 1 we  can  use  (2.16)  to  obtain 


P { max  = exp( e" ^ 

1 < j < n ^ 


where  a = log  p"^  and  = a"^  log(n  s^(  l-p)/s.') , Note  that  this  is 
consistent  with  (4.3).  ^ 


6.  Numerical  Results 

A simulation  of  the  I queue  was  performed  to  help  assess  the 

effectiveness  of  the  estimation  methods  proposed  in  Section  3.  0“*^  goal 

is  to  estimate  the  d.f.  of  ^ n where  is  either  the 

maximum  waiting  time^  virtual  waiting  time^  or  queue- length  in  the  kth 

cycle  and  n is  the  number  of  cycles  simulated.  These  d. f. 's  will  all 

be  estimated  by  A( (x-b^)/a^) . Thus  our  principal  task  is  to  estimate 

the  appropriate  a^  and  b^  for  the  three  processes  mentioned  abovej 

theoretical  values  of  a and  b are  available  from  Section  4. 

n n 
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The  random  number  generator  used  was  subroutine  GGU5  available 

in  the  IMSL  package.  This  generator  is  the  congruential  generator  developed 

by  LEARMONTH  and  LEWIS  (19^73).  Regression  lines  were  fitted  to  the  plots 

of  log  E^(x)  versus  x beginning  at  x = 0.  For  the  ^ 

(V^:t  > 0}  processes  observations  were  grouped  in  I50  cells  of  length 

0.02  [0.05]  for  p = 0.5  [0.9].  Tables  L-9  contain  the  results  of  the 

simulation  for  estimating  a^  and  b^.  Two  values  of  p (O.5  and  O.9) 

were  used  along  with  several  combinations  of  run  lengths  (number  of 

cycles  simulated)  and  number  of  replications.  The  entries  contained  in 

the  tables  are  the  sample  means  of  the  various  estimates  over  the  number 

of  replications  and  the  half-length  of  a symmetric  90^  confidence 

interval  about  the  sample  mean.  For  example,  in  Table  4 take  p = 0«5, 

500  cycles,  and  100  replications.  Then  for  k = 20  in  the  Feigin-Weissman 

♦ 

procedure  0. 1899  is  the  sample  mean  of  100  estimates,  ^qq,  ^500 
[.1899  - .0077,  .1899  + .0077]  is  the  corresponding  90^:  confidence 
interval  based  on  the  100  estimates.  In  this  case  the  true  value  of 


a is  0.2. 
n 

Here  are  some  general  observations  based  on  Tables  14--9. 

AUMVUE  estimacors  a and  b do  in  general  have  smaller  bias  than 

n n 

/s  ^ 

the  AMLE  estimators  a and  b . however,  the  mean  square  errors  (mSE)  , 
that  is  bias  squared  plus  variance,  are  comparable.  In  general,  the 
bias  decreases  with  increasing  k.  The  regression  estimator  in  general 
yields  a smaller  bias  than  the  other  estimators  for  the  ^ 

(Vj,:t  > 0)  processes.  In  those  cases  where  the  bias  is  larger  the  MSE 
is  usually  quite  a bit  smaller.  For  the  (Q^:t  > 0)  process  the  smallest 
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TABLE  L 


Estimates  of  for  > 0)  in  the  M/M/ 1 queue  with  4 


# cycles/ 

# repL. 

P 

true 

value 

regression 

a 

n 

k = 

10 

k = 

20 

k = 30 

a 

n 

♦ 

a 

n 

/% 

a 

n 

♦ 

a 

n 

a 

n 

♦ 

a 

a 

250/200 

0-5 

B 

.20L6 

. 1678 

.1864 

.1807 

.1903 

.1799 

.1861 

m 

.0059 

.0070 

.0078 

.0052 

.0055 

.0041 

.00ii3 

500/ 100 

0.5 

.2008 

. 1704 

.1893 

. 1804 

.1899 

.1818 

.1881 

m 

.0065 

.0090 

.0100 

.0073 

.0077 

.0063 

.0065 

1000/50 

0.5 

.1972 

.1721 

.1912 

.1832 

.1928 

.1793 

.1855 

B 

.0068 

.0109 

.0121 

.0068 

.0093 

.0081 

.0083 

1000/80 

0.9 

1.0 

.&Jh6 

.8172 

.9080 

.8313 

•8751 

.8207 

.8490 

.0398 

.0508 

. C^64 

.0391 

.0411 

.0319 

.0330 

2000/ UO 

0.9 

1.0 

.8739 

.8515 

.9461 

.8667 

.9123 

.8636 

.893^^ 

.0486 

.0821 

.0912 

.0607 

.0639 

.0500 

.0517 

I4.OOO/20 

0.9 

1.0 

.8865 

.8823 

.9803 

•9312 

.9802 

.8984 

.9293 

.0591 

.1298 

.1442 

.0895 

.0942 

.0773 

.0799 

26 


TABLE  5 


Estimates  of  for  > 0}  fn  the  M/Vl  queue  with  = 10 


1 

1 

1 

k = 

10 

k = 

20 

k = 30 

# cycles' 

# repl. 

P 

true 

value 

regres- 

sion 

b 

n 

/s 

b 

n 

* 

b 

n 

b 

n 

* 

b 

n 

/s 

b 

n 

♦ 

b 

n 

250/200 

0.5 

0.8270 

VO  00 
ovoo 

•-> 

00  0 

0 d 

0.77^ 

0.0178 

0.8061 

0.0190 

0.7883 

0.0162 

0.8120 

0.0168 

0.781+8 

0.011+6 

0.3028 

0.0150 

500/  100 

0.5 

0.9657 

0.9761+ 

0.0239 

0.9068 

0.021+7 

0.91+08 

0.0263 

0.9189 

0.0236 

0.91+25 

0.021+5 

0.9201 

0.0221 

0.9382 

0.0227 

1000/50 

0.5 

1. 101+5 

1. 1016 
0.0280 

1.01+21+ 

0.0285 

1.0767 

0.0501+ 

i.055i^ 

0.0277 

1.0795 

0.0287 

1.01+56 

0.0281+ 

1.0635 

0.0292 

1000/80 

1 

1 

0.9 

i+.  1+998 

1+.  1+067 
0.11+98 

1+ . 1663 
0. 11+1+6 

1+.3292  ' 

0.151+1 

i+.  1605 
0. 1512 

4.2695 

0.1361 

1+.1511 

0.1201+ 

1+.2131 

0.1255 

2000/1+0 

1 

5.1950 

1 

5.0102 

0.2170 

!+.8i+7ii- 

0.21+60 

5.0172 

0.2611+ 

i+.  81+61+ 
0.2182 

I+.960O 

0.2259 

1+.8520 

0.1967 

4.91S3 

0.2015 

1+000/20 

1 

5.8861 

5.6715 

0.2812 

5.5718 

0.3362 

5.7J^77 

0.3598 

5.6337 

0.3089 

5.7558 

0.3202 

5.5565 

0.291+3 

5 .61+65 
0.3019 

TABLE  7 


Estimates  for  'b  for  (V  :t  > 0}  in  the  M/M/1  queue  with  ^ = 10 


# cycles/ 

# repl. 

1 

regres- 

sion 

b 

n 

k = 

10 

1 k = 

1 

20 

k = 

50 

P 

true 

value 

b 

n 

b* 

n 

b^ 

b 

n 

b 

n 

* 

b 

n 

250/200 

0,5 

0.9657 

I 

0.9797 

0.0180 

0.9140 

0.0178 

0.9477 

0.0191 

0.9257 

0.0162 

0.9492 

0.0168 

0.9295 

0.0150 

0.9475 

0,0154 

“SOO/  100 

0.5 

1. 1045 

1. 1103 
0.0240 

1.0428 

0.0243 

1.0762 

0.0259 

1.0550 

0.0226 

1.0785 

0.0234 

1.0606 

0.0219 

1.0788 

0.0224 

1000/50 

0.5 

1.2L29 

1.2584 

0.0304 

1.18?7 

0.0555 

'1.2171 

0.0359 

1. 1889 
0.0305 

1.2125 

0.0317 

1.1975 

0.0283 

1.2157 

0.0291 

1000/80 

0.9 

4.6052 

4.4714 
0. 1474 

4.2752 

0.1455 

1 

4.4383 

0.1555 

4.2645 

0.1519 

4.3755 
0. 1368 

4.2278 
0. 1206 

4.3094 

0.1237 

2000/14-0 

0.9 

5.2983 

5.0695 

0.2117 

4.9571 

0.2450 

5.1^1 

0.2602 

4.9495 

0.2210 

5.0628 

0.2290 

4.9512 

0.1959 

5.0380 

0.2006 

I4.OOO/20 

0.9 

5.9915 

5.6715 

0.2812 

5.6769 

0.3420 

5.8550 

0.5654 

5 .7259 
0.5126 

5.8447 

0.3242 

5.6716 

0.2977 

5.7616 

0.3055 

1 


TABLE  8 

Estimates  of  for  > 0)  in  the  M/M/l  queue  with  = 10 


1 

# cycles/ 

# repl. 

i 

k = 

10 

k = 

20 

k = 

30 

P 

true 

value 

regres- 

sion 

a 

n 

/S 

a 

n 

* 

a 

n 

/N 

a 

n 

1 

* 

a 

n 

a 

n 

* 

a 

n 

250/200 

0.5 

i.i+i+2r 

1.5^01+ 

0.0317 

1 

1.251+1+ 

0.0566 

1.5958 

0.0629 

1 

1.55i^9 

0.01+77 

1.4262 

0.0502 

1.15?7 

0.0404 

1.1925 

0.0418 

500/100 

0.5 

1.1+1+27 

1.5538 

0.0567 

1.1919 

0.071+8 

1.321+1+ 

0.0851 

1.3600 

0.0668 

1.4515 

0.0703 

1.2116 

0.0613 

1.2534 

0.0635 

1000/50 

0.5 

1.1+1+27 

1.3297 

0.0395 

1.2020 
0. 1120 

1.3355 

0.121+1+ 

1.361+0 

0.0903 

1.4558 

0.0950 

1.5126 

0.0923 

1.3579 

0.0949 

1000/80 

0.9 

9.1+912 

7.9711+ 

0.1+069 

7.6221+ 

0.5067 

8 . 1+691+ 
0.5630 

7.7721+ 

0.3883 

8.1815 

0.4088 

7.8412 

0.3^79 

8.1116 

0.3392 

2000/1+0 

0.9 

9.1+912 

7.9621 

0.5117 

8.2575 

0.81+11 

9. 1750 
0.93^5 

8.01+62 

0.61+18 

8.4697 

0.6756 

8 . 0466 
0.5252 

8.3241 

0.5435 

1+000/20 

0.9 

9.1+912 

8.21+06 

0.6612 

8.7600 

1.2651 

9.7535 

1.1^057 

8.1+350 

0.8397 

8.8789 

0.8839 

8.6533 

0.7252 

8.9517 

0.71+81 

1 

i 


TABLE  9 


Estimates  of  for  (Q^:t  > 0)  in  the  IV^m/1  queue  with 


# cycles/ 


regres 

Sion 


0.5  6.9658  6.9196  7.0985  7.5^6  7.2691  7.41^68  6.8159  6.9511 

0.1222  0.1315  O.lJ+17  0.1278  0.1358  0.1202  0.121;1 

7.9658  7.8UI  7.97*^  8.2122  8.201+1  8.3825  7.8811  8.0021 

0.1603  0.1769  0.1900  0.1890  0.1971+  0.1727  0.1786 

8.9658  8.7129  8.9077  9.11+73  9.1611  9.3^50  9.061+6  9.1958 

0.1911  0.21+56  0.2657  0.2558  0.2671  0.21+57  O.25L5 


1+0.7508  1+2.2701+  1+0.6961  1+1.7151+  1+0.7690  lil.5525 
IA567  1.5525  1.2861  1.3355  1.1925  1.221+5 

1+7.5382  1+9.181+1+  1+7.0291  1+8.081+1  1+6.9179  1+7.7221 
2.5152  2.671+6  2.2171+  2.2985  2.0136  2.061+6 


3.1^052  3.61+25  2.9162  3.0220  2.71+75  2. 3183 


0.9 

1+3.7087 

1 

50.2875 

1 

56.8665 

bias  when  estimating  a was  obtained  from  a with  k = 10  or  20. 

n n 

Again^  however  the  MSE  of  the  regression  estimator  was  generally  smallest. 
When  estimating  b^  for  > 0}  the  smallest  bias  was  obtained  by 

four  estimators  in  different  cases.  Percentage-wise  the  differences  in 
bias  were  small  and  again  the  MSE  of  the  regression  estimator  was  often 
the  lowest. 

* * 

To  summarize  these  results,  we  recommend  using  the  a and  b 

* n n 

estimators  rather  than  a and  b . The  value  of  k should  be  as  large 

n n 

as  possible,  realizing  of  course  that  the  amount  of  computer  time  required 

will  increase  with  k.  The  regression  estimators  a^  and  b^  performed 

very  well.  They  achieved  the  smallest  bias  in  about  one- third  of  the  cases 

and  had  the  smallest  MSE  in  virtually  all  cases.  In  practice,  we  suggest 

that  a plot  of  log  E^(x)  versus  x be  made  before  the  regression  line 

is  fitted.  This  will  enable  the  simulator  to  eyeball  a straightline  fit 

and  the  x-value  at  which  to  begin  fitting. 

Our  ultimate  objective  is  to  estimate  the  extreme  value  distribution 

P{  max  < x)  for  the  three  queueing  processes.  Each  of  these 

distributions  will  be  approximated  by  A((x-b^)/a^)  and  then  the  constants 

a and  b estimated  by  one  of  the  methods  discussed  above.  The  approxima 
n n 

tion  of  P(  max  ^ ^ A((x-b  )/a  ) is  a large  sample  one  valid 

l<k<n^  "" 

for  large  n from  the  limit  theorem  (2.1).  As  the  simulation  results  for 
the  waiting  time  and  virtual  waiting  time  processes  were  quite  similar,  we 
only  present  the  figures  for  the  waiting  time  and  queue  length  processes. 

For  the  Feigin-Weissman  (F-W)  method  only  the  k = 30  case  is  presented. 
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TABLE  10 


Estimates  of  P{  max  < x]  for  (W  :n  > 0} 

1 < k < n ^ 

in  the  M/M/1  queue,  p = 0.5,  u = 10 


# cycles/ 

method  | § repl. 


regression  250/20O 

30  2^0l200 


A((x-bJ/a^)/x 


regression  500/ 1-00 
F-W,k  = 30  500/100 


1000/50 

1000/50 


.25/ .76 17 

.50/. 9003 

.75/1.0762 

.90/1.2771 

.99/1.7471 

.2876 

.4894 

.6994 

.8436 

.9721 

.0261 

.0289 

.0246 

.0168 

.0054 

.5358 

.5578 

.7666 

.8983 

.9867 

.025k 

.0256 

.0192 

.0112 

.0023 

A((x.bn)/ajj)/ 

X 

.25/. 9005 

.50/1.0390 

.75/1.2148 

.90/1.4157 

.99/1.8857 

.2870 

.4962 

.7148 

.8647 

.9788 

.0552 

.0373 

.0354 

.0199 

.0054 

.3448 

.5586 

.7637 

.8946 

.9858 

.0387 

.0380 

.0280 

.0161 

.0032 

A((x-b  )/aJ/ 

X 

.25/1.0390 

.50/1.1776 

.75/1.3535 

.90/1.5544 

.99/2.0243 

.2921 

.5147 

.7382 

.8832 

.9845 

.0422 

.0459 

.0570 

.0226 

.0047 

.5615 

.5820 

.7842 

.9071 

.9882 

.0512 

.0501 

.0361 

.0202 

.0038 

TABLE  11 

Estimates  of  P{  max  < x)  for  (W  :n  > 0) 

1 < k < n ^ " 


in  the  M/M/ 1 queue^  p = 0.9,  4 = 10 


method 

# cycles/ 

# repl. 

.25/4.1752 

A((x-b^)/a^)/x 

.50/4.8663  .75/5.71^57 

.90/6.7502 

.99/9.1000 

regression 

1000/30 

.5i^95 

.0475 

.5654  .7674 

.0484  .0380 

.8959 

.0248 

.9837 

.0068 

F-W^k  = 30 

1000/80 

.4003 

.0437 

.6241  .8141 

.0423  .0308 

.9225 

.0175 

.9906 

.0052 

A((x-bj/aj/x 

.25/4.8663 

.50/5.5595  .75/6.489 

.90/7.1+1+53 

.99/9.7951 

regression 

2000/1+0 

.5732 

.0627 

.591+8  .7901 

.0628  .0520 

.9046 

.0579 

.981+2 

.0120 

F-W^k  = 30 

2000/1+0 

.4074 

.0701 

.6037  .7719 

.0716  .0613 

.8788 

.0461 

.9702 

.0194 

i 

A((x-bj/aj/x 

.25/5.5595 

.50/6.2526  .75/7.1320 

.90/8.1365 

.99/10.4862 

regression 

4000/20 

.5691 

.0851 

.5956  .7923 

.0839  .0691 

.9094 

.0449 

.9877 

.0099 

F-W,k  = 30 

4000/20 

.3881 

.0942 

.5928  .7798 

.0924  .0720 

.8979 

.0450 

.981+2 

.0105 

TABLE  12 

Estimates  of  P(  max  ^ (Q  :t  > 0) 

l<k<n^  ^ ~ 

in  the  M/M/ 1 queue,  p = 0.5,  ^ = 10 


method 

flf  cycles/ 

# repl. 

.57^7 

P{  max  < 

1 < k < n ^ 

.6128/8  .7831/9 

x]/x 

.9408/11 

.9924/14 

regression 

250/200 

.4289 

.0276 

.6517  .7785 

.0261  .0208 

.9272 

.0105 

.9868 

.0051 

F-W,k  = 50 

250/200 

.hk3l 

.0514 

.6529  .7975 

.0281  .0218 

.9358 

.0106 

.9884 

.0050 

P(  max  M,"^  < 

1 < k < n ^ ~ 

x]/x 

.3755/8 

.6152/9  .7852/10 

.9408/ 12 

.9924/15 

regression 

500/100 

.4451 

.0569 

.6494  .7944 

.0545  .0265 

.9364 

.0121 

.9897 

.0029 

F-W,k  = 50 

500/ 100 

.4251 

.0447 

.6505  .7781 

.0404  .0515 

.9264 

.0156 

.9860 

.0043 

P(  max  m7  < 

1 < k < n ^ 

x)/x 

.6155/10  .7855/11 

.9408/15 

.9924/  16 

regression 

1 

1000/50 

.4621 

.0462 

.6707  .8156 

.0419  .0511 

.9467 

.0127 

.9924 

.0025 

F-W,k  = 50 

1000/50 

.3325 

.0582 

.5891  .7427 

.0573  .0480 

.9065 

.0261 

.9801 

.0080 

TABLE  15 


Estimates  of  P(  max  ^ ^ (Q  ;t  > 0) 

1 < k < n 

in  the  M/M/1  queue,  p = O.9,  n = 10 


method 

§ cycles/ 

# repl. 

.2599/ LI 

P(  max  m7  < 

1 < k < n ^ 

.5&72/L8 

x)/x 

.9088/66 

.9906/88 

regression 

1000/80 

.3907 

.0506 

.6156  .79?/ 

.0L9I  .0388 

.9080 

.0262 

.9836 

.0088 

\ 

F-W,k  = 30 

1000/80 

.399^^ 

.0L36 

.6336  .8101 

.OL33  .0351 

.9215 

.0189 

.9893 

.0038 

f 

P[  max  m7  < 

1 < k < n ^ 

x}/x 

.2800/L8 

.5O7L/5L  .7693/63 

.903L/72 

.9900/94 

regression 

2000/ lj.0 

>396 

.0700 

.6321  .8207 

.0657  .0531 

.9126 

.04.10 

.9830 

.0144 

F-W,k  =30 

2000/ I4.O 

.i^379 

.0697 

.6233  .811L 

.0661  .0505 

.9093 

.0344 

.9849 

.0096 

i 

P(  max  M,^  < 

1 < k < n ^ 

x}/x 

. 2575/5  iv 

.5231/61  .7568/69 

.9074/79 

.9905/101 

regression 

LOOO/20 

.i^033 

.0905 

.6318  .8027 

.0902  .0790 

.9116 

.0537 

.9853 

.0133 

F-W,k  = 30 

1^.000/20 

.5712 

.0968 

.5860  .7696 

.0906  .0691 

.8978 

.0421 

.9842 

.0104 

the  results  in  general  being  better  than  for  smaller  k.  As  an  example 
of  the  entries  in  Tables  10-15,  consider  in  Table  10  the  regression  estimate 
for  500  cycles  and  100  replications  when  A((x-b^)/a^)  = 0.75  with 
X = 1.2148.  The  mean  of  the  100  estimates  in  this  case  is  0.7148  and  the  90% 
confidence  interval  has  half-length  0.0554-.  For  the  queue  length  process 
the  exact  value  of  P{  max  ^ I’®  calculated  from  the  results 

in  Example  (5.4-),  and  these  values  are  indicated  in  Tables  12  and  I5.  A 
comparison  of  the  bias  and  MSE  for  the  regression  and  F-W  estimates  in 
Tables  10- I5  shows  them  to  be  quiet  comparable.  For  both  estimators  the 
bias  is  in  general  smaller  for  the  larger  percentiles.  For  the  smallest 
percentiles  both  estimators  over  estimate  the  quantity  in  question. 

Recommendations 

The  first  step  in  estimating  extreme  values  for  regenerative  simula- 
tions is  to  plot  the  function  log  E^(x)  versus  x.  This  enables  the 
simulator  to  ascertain  whether  the  tail  of  the  distribution  of  the  maximum 
value  within  a cycle  is  exponential.  Based  on  the  examples  we  can  calculate, 
we  would  expect  this  tail  to  be  exponential  for  regenerative  processes 
with  finite  mean  cycle  length.  If  the  plot  of  log  E^(x)  is  roughly 
linear,  we  would  then  recommend  estimating  a and  b by  both  the 

regression  and  Feigin-Weissman  AUMVUE  methods.  For  the  F-W  method  k 
should  be  at  least  equal  to  50,  and  perhaps  even  larger  if  increased 
computer  time  is  not  a big  problem.  We  would  expect  the  two  methods  to 
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give  farily  comparable  results.  Using  both  methods  gives  the  simulator 

a rough  check  against  inadvertant  programming  errors.  Based  on  our 

simulations  of  the  M/M/ 1 queue  we  would  select  the  regression  method 

when  the  program  is  debugged^  if  forced  to  choose  just  one  method.  Once 

estimates  for  a and  b are  available,  the  distribution  of  extreme 
n n > 

values  for  n cycles  (P{  max  M^  < x))  can  be  approximated  by 

1 < k < n ^ 

A((x-b^)/an)  where  the  appropriate  estimates  of  a^  and  b^  are  used. 
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"Regenerative  Simulation  for  Extreme  Values"^  Technical  Report  No.  45. 

or  - 

^Let  {X^:t  >0}  denote  the  regenerative  process  being  simulated  and 

assume  that  converges  weakly  (in  distribution)  to  a limit  random 

variable  X.  Our  conern  itr~thira~-pap«r  is  in  estimating  the  extreme 

values  of  the  process  {X^:t  >^0].  Suppose  we  are  interested  in  the 

largest  value  attained  in  the  interval  (O, t);  X^  = sup<X^:0  < s < t>. 

Examples  of  this  are  the  maximum  queue  lengths  or  waiting  times  in  a 

queueing  system.  As  t increases  so  will  X^^  without  bound  if  the 

state  space  of  {Xj.:t  > 0)  is  unbounded.  This  report  developes  two 

* 

methods  for  estimating  the  distribution  of  X^.  When  the  regenerative 

process  is  either  the  GI/g/1  queue  or  a birth-death  process  theoretical 

* 

results  are  available  for  the  distribution  of  X^.  The  waiting  time^ 

queue  lengthy  and  virtual  waiting  time  for  an  m/m/ 1 queue  were  simulated. 

* 

The  two  methods  for  estimating  the  distribution  of  X^  were  employed 
and  the  simulation  results  compared  with  the  theoretical  results. 
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